function dXdt = hrfnhr(t, X, eps_val, params)
% 计算 HR-FN-HR 神经网络模型中各状态变量的导数
%
% 状态变量顺序：X = [x1; x2; x3; x4; x5; x6; phi]
% 其中：
%   x1  - HR 神经元1的膜电位
%   x2  - HR 神经元1的恢复变量
%   x3  - FN 神经元的膜电位
%   x4  - FN 神经元的恢复变量
%   x5  - HR 神经元2的膜电位
%   x6  - HR 神经元2的恢复变量
%   phi - memristor 状态变量
%
% 微分方程如下（参数均来自 params 结构体，eps_val 为当前 ε）：
%
%   dx1/dt = x2 - a1*x1^3 + b1*x1^2 + i1 + m12*tanh(phi)*(x3 - x1)
%   dx2/dt = c1 - d1*x1^2 - x2
%   dx3/dt = x3 - b2*x3^3 - x4 + i2 + m21*(x1 - x3) + m23*tanh(phi)*(x5 - x3)
%   dx4/dt = eps_val*(a2 + x3 - c2*x4)
%   dx5/dt = x6 - a3*x5^3 + b3*x5^2 + i3 + m32*(x3 - x5)
%   dx6/dt = c3 - d3*x5^2 - x6
%   dphi/dt = -0.2*phi^3 + x1*x3*x5 + phi
%
% 输入:
%   t       : 时间（本模型中 t 未直接使用，但必须保留以符合 ODE 求解器接口）
%   X       : 状态向量 [x1; x2; x3; x4; x5; x6; phi]
%   eps_val : 当前 FN 神经元参数 ε 的值
%   params  : 参数结构体，包含以下字段：
%               a1, b1, c1, d1, i1         % HR 神经元1参数
%               a2, b2, c2, i2              % FN 神经元参数
%               a3, b3, c3, d3, i3          % HR 神经元2参数
%               m12, m21, m23, m32          % 耦合系数
%
% 输出:
%   dXdt    : 各状态变量的导数，排列为 [dx1/dt; dx2/dt; dx3/dt; dx4/dt; dx5/dt; dx6/dt; dphi/dt]

    % 虽然 t 在本模型中未直接使用，但为了满足 ode45 接口保留 t 参数
    % 将状态变量展开，便于理解各变量
    x1  = X(1);
    x2  = X(2);
    x3  = X(3);
    x4  = X(4);
    x5  = X(5);
    x6  = X(6);
    phi = X(7);
    
    % 分别计算各变量的导数
    dx1_dt = x2 - params.a1*x1^3 + params.b1*x1^2 + params.i1 + params.m12*tanh(phi)*(x3 - x1);
    dx2_dt = params.c1 - params.d1*x1^2 - x2;
    dx3_dt = x3 - params.b2*x3^3 - x4 + params.i2 + params.m21*(x1 - x3) + params.m23*tanh(phi)*(x5 - x3);
    dx4_dt = eps_val*(params.a2 + x3 - params.c2*x4);
    dx5_dt = x6 - params.a3*x5^3 + params.b3*x5^2 + params.i3 + params.m32*(x3 - x5);
    dx6_dt = params.c3 - params.d3*x5^2 - x6;
    dphi_dt = -0.2*phi^3 + x1*x3*x5 + phi;
    
    % 将所有导数组合成列向量返回
    dXdt = [dx1_dt; dx2_dt; dx3_dt; dx4_dt; dx5_dt; dx6_dt; dphi_dt];
end
